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REGULARIZATION  OF  THE  CHAPMAN-ENSKOG  EXPANSION  AND 
ITS  DESCRIPTION  OF  SHOCK  STRUCTURE* 


KUN  XU+ 

Abstract.  In  the  continuum  transition  flow  regime,  we  propose  to  truncate  the  Chapman- Enskog 
expansion  of  the  Boltzmann  equation  to  the  Navier-Stokes  order  only  without  going  to  the  Burnett  or  super 
Burnett  orders.  However,  the  local  particle  collision  time  has  to  be  generalized  to  depend  not  only  on  the 
local  macroscopic  flow  variables,  but  also  their  gradients  in  the  rarefied  gas  regime.  Based  on  the  gas-kinetic 
BGK  model,  the  relation  between  the  conventional  collision  time  and  the  general  one  is  obtained.  More 
speciflcially,  a  generalized  constitutive  relation  for  stress  and  heat  flux  is  proposed.  This  new  model  is 
applied  to  the  study  of  argon  gas  shock  structure.  There  is  good  agreement  between  the  predicted  shock 
structure  and  experimental  results  for  a  wide  range  of  Mach  numbers. 

Key  words.  Navier-Stokes  equations,  Chapman-Enskog  expansion,  continuum  transition  flow,  shock 
structure 

Subject  classification.  Applied  Numerical  Mathematics 

1.  Introduction.  It  is  well  recognized  that  the  Navier-Stokes  equations  of  the  classical  hydrodynamics 
are  incapable  of  accurately  describing  shock  wave  phenomena  and  also  for  the  flow  phenomena  in  the  rarefied 
regime.  In  order  to  improve  the  Navier-Stokes  solutions,  much  effort  has  been  paid  on  the  construction  of 
higher-order  hydrodynamic  equations  based  on  the  Chapman-Enskog  expansion.  But  the  Burnett  and  super 
Burnett  equations  give  unstable  shock  structures  in  high  Mach  number  cases.  Eor  example,  no  shock  structure 
can  be  obtained  for  the  Burnett  equations  after  a  critical  Mach  number  =  2.69  [9].  Even  though  the 
argumented  Burnett  of  Zhong  et  al.  and  BGK-Burnett  equations  of  Agarwal  et  al.  can  significantly  improve 
the  Navier-Stokes  solutions  in  the  continuum  transition  regime  [7,  14,  1],  it  is  unclear  that  the  stable  shock 
structures  of  these  schemes  are  coming  from  the  complicated  numerical  dissipations,  such  as  the  use  of 
Steger  Warming  flux  splitting  scheme  for  the  inviscid  part  of  the  equations  [8] ,  or  the  selected  higher-order 
terms.  As  analyzed  in  [10],  the  failure  of  the  Burnett  equations  for  the  shock  structure  calculation  is  not  too 
surprising  because  the  applicability  of  the  Chapman-Enskog  expansion  itself  is  valid  to  the  small  Knudsen 
numbers.  The  possible  generation  of  spurious  solutions  from  the  higher-order  terms  in  the  Chapman-Enskog 
expansion  is  another  point  for  criticism  [5]. 

This  work  is  motivated  originally  by  extending  the  gas-kinetic  BGK  Navier-Stokes  solver  to  the  con¬ 
tinuum  transition  regime  [12].  The  direct  adoption  of  the  Chapman-Enskog  expansion  with  the  terms  pro¬ 
portional  to  Knudsen  number  and  in  the  gas  distribution  function  encounters  great  difficulty  in  the 
shock  structure  calculations.  The  critical  Mach  number  for  the  shock  structure  based  on  the  BGK-Burnett 
expansion  is  found  to  be  around  =  4.5,  and  the  number  becomes  even  smaller,  i.e.,  =  2.0,  with  the 

inclusion  of  super  Burnett  term  [13].  Our  numerical  experiments  show  clearly  that  the  successive  Chapman- 
Enskog  expansion  without  selectively  choosing  higher  order  terms  give  divergent  results  as  the  Knudsen 
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support  was  provided  by  Hong  Kong  Research  Grant  Council  through  RGC  HKUST6132/00P. 

tMathematics  Department,  Hong  Kong  University  of  Science  and  Technology,  Clear  Water  Bay,  Kowloon,  Hong  Kong 
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number  increases.  However,  up  to  the  Navier-Stokes  order,  there  is  not  any  limitation  on  the  Mach  number 
for  the  existence  of  the  stable  shock  structure.  This  observation  is  consistent  with  the  theoretical  analysis 
in  [9].  Therefore,  it  may  be  possible  to  truncate  the  Chapman- Enskog  expansion  to  the  Navier-Stokes  order 
only  and  include  the  possible  non-equilibrium  effect  on  the  modification  of  the  viscosity  and  heat  conduction 
coefficients,  the  so-called  constitutive  relations.  Traditionally,  the  particle  collision  time  r  is  regarded  as  a 
function  of  macroscopic  variables.  For  example,  based  on  the  BGK  model  [3],  we  have  the  collision  time 
r  =  /i/p,  where  /i  is  the  dynamical  viscosity  coefficient,  such  as  the  Sutherland’s  law,  and  p  is  the  pressure. 
All  those  viscosity  coefficients  are  basically  obtained  either  experimentally  or  theoretically  in  the  continuum 
flow  regime  [6].  There  is  no  reason  to  guarantee  that  this  relation  is  still  applicable  for  the  rarefied  gas.  In 
this  paper,  we  are  going  to  derive  a  general  particle  collision  time  r*,  which  is  applicable  in  both  continuum 
and  continuum  transition  regime.  This  derivation  is  based  on  the  closure  of  the  Chapman- Enskog  expansion 
on  the  Navier-Stokes  order  and  the  BGK  equation. 

2.  Closure  of  the  Chapman-Enskog  Expansion.  The  BGK  model  in  the  x-direction  can  be  written 
as  [5] 

(2.1)  ft+uf,  =  ^, 

T 

where  /  is  the  gas  distribution  function  and  g  is  the  equilibrium  state  approached  by  /.  Both  /  and  g 
are  functions  of  space  x,  time  t,  particle  velocities  u,  and  internal  variable  The  particle  collision  time  r 
determines  the  viscosity  and  heat  conduction  coefficients,  i.e.,  p  =  rp.  The  equilibrium  state  is  a  Maxwellian 
distribution, 

TT 

where  p  is  the  density,  U  is  the  macroscopic  velocity  in  the  x  direction,  and  A  is  related  to  the  gas  temperature 
to/2A:T.  For  a  monatomic  gas,  and  ^2  represent  the  particle  velocities  in  the  y  and  2:  directions.  The 
relation  between  mass  p,  momentum  pU,  and  energy  pE  densities  with  the  distribution  function  /  is 

I  p[7  I  =  / ‘tpfdudiidi‘2, 

\pe] 

where  -i/j  has  the  components 

Since  mass,  momentum  and  energy  are  conserved  during  particle  collisions,  /  and  g  should  satisfy  the 
compatibility  condition 

(2.2)  I (g  -  f)xPadud^^d^2  =  0,  a  =  1,2,3, 

at  any  point  in  space  and  time. 

It  is  well  known  that  the  Euler,  the  Navier-Stokes,  and  the  Burnett,  etc.  equations  can  be  derived 
from  the  above  BGK  model  using  the  Chapman-Enskog  expansion  [6].  The  successive  expansion  of  the 
Chapman-Enskog  expansion  gives 

f  —  9  T  '^9x^  T  {dtt  T  T  ^  9xx^  {dttt  T  ^'^9xtt  T  gxxt  T  ^  9xxx^  T  ■■■ 
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which  corresponds  to  the  Euler  (r°),  the  Navier-Stokes  (r),  the  Burnett  (r^),  and  the  super  Burnett  (r^)  ... 
orders.  With  the  definition  D  =  d/dt  +  ud/dx,  we  can  write  the  above  equation  as 


n— 1 

In  the  continuum  transition  regime,  the  Navier-Stokes  equations  are  expected  to  be  inaccurate  and  the 
expansions  beyond  the  Navier-Stokes  order  have  only  achieved  limited  success.  As  shown  by  Uribe  et  al. 
[10],  Bobylev’s  instability  analysis  basically  provides  a  range  of  Knudsen  numbers  for  which  the  Burnett 
order  is  valid  [4]. 

In  order  to  increase  the  validity  of  the  gas  kinetic  approach  in  the  continuum  transition  regime,  we  have 
to  regularize  the  Chapman-Enskog  expansion.  The  main  idea  of  this  paper  is  to  close  the  Chapman-Enskog 
expansion  up  to  the  Navier-Stokes  order  only  without  going  to  Burnett  or  super  Burnett  orders.  But,  instead 
of  keeping  the  original  particle  collision  time  r,  we  have  to  construct  a  general  one.  In  other  words,  we  expand 
the  gas  distribution  function  as 


(2.3) 


f  =  9-  T*{9t  +ugx), 


and  r*  is  obtained  to  have  the  BGK  equation  to  be  satisfied. 


(2.4) 


/  =  9  -r{ft+uf^). 


When  the  spatial  and  temporal  derivatives  of  the  particle  collision  times  are  ignored,  from  the  above  two 
equations  (2.3)  and  (2.4),  we  can  get  the  relation  between  the  original  particle  collision  time  r  and  the  new 
one  r* , 

(2.5) 


l+rD^glDg- 

Therefore,  the  local  particle  collision  time  depends  not  only  on  the  macroscopic  variables  through  r  =  /i/p, 
but  also  the  ratio  between  the  Burnett  order  D^g  and  the  Navier-Stokes  order  Dg.  In  the  above  equation, 
r*  depends  on  the  particle  velocities,  which  may  introduce  great  complexity  in  using  its  solution.  In  order 
to  remove  the  particle  velocity  dependence  in  r*,  we  suggest  to  take  a  moment  on  D‘^g/Dg  first,  such  as 


(2.6) 


_ 
Dg 


)  =  J  '^{u)D^gdud^id^2/  j  '^{u)Dgdud^id^2- 


Here  we  propose  to  use  'l/(u)  =  {u  —  in  the  above  integration,  where  U  is  the  local  macroscopic  velocity. 
Other  choices  may  be  possible.  But,  due  to  the  fact  that  both  moments  of  Dg  and  D‘^g  on  (l,u,  (l/2)(u^  + 
+  ^D)  vanish,  the  above  choice  becomes  the  only  one  which  mimics  ‘dissipative’  energy  in  some  sense. 
In  the  expressions  of  D‘^g  and  Dg,  there  exist  temporal  and  spatial  derivatives  of  a  Maxwellian.  The  local 
spatial  derivatives  can  always  be  constructed  from  the  interpolated  macroscopic  flow  variables,  such  as  the 
gradients  of  mass,  momentum,  and  energy.  Eor  the  temporal  derivatives,  they  have  to  be  evaluated  based 
on  the  compatibility  conditions,  such  as  f  D"^ g'tpadud^id^^  =  0  and  f  Dg'tpadud^id^^  =  0  of  the  Chapman- 
Enskog  expansion.  The  detailed  numerical  procedure  is  given  in  [13].  In  summary,  based  on  the  BGK  model 
and  the  closure  of  the  Chapman-Enskog  expansion  on  the  Navier-Stokes  order,  we  derive  a  new  local  particle 
collision  time  r* ,  such  that 


(2.7) 


n  =  r/  (1  -b  T{D‘^g/Dg))  . 
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Based  on  the  above  r*,  the  viscosity  and  heat  conduction  coefficients  will  depend  on  both  the  macroscopic 
variables  and  their  slopes.  In  the  continuum  regime,  since  the  higher-order  dissipation  should  have  less  effect 
than  the  lower  order  one,  (D^g/Dg)  will  theoretically  go  to  zero.  This  is  verified  in  the  following  shock 
structure  calculation. 

In  recent  years,  an  accurate  gas-kinetic  BGK  Navier-Stokes  solver  (BGK-NS)  has  been  developed  for  the 
viscous  solution  in  the  continuum  regime  by  the  current  author  and  co-workers  [12] .  In  the  following  argon  gas 
shock  structure  calculations,  we  are  going  to  use  the  above  BGK-NS  method,  but  with  the  implementation  of 
the  new  particle  collision  time  r*.  For  a  monatomic  gas  modeled  by  point  centers  of  force,  the  kinetic  theory 
leads  to  a  viscosity  /x  proportional  to  and  the  Prandtl  number  Pr  =  i^Cp/n  is  a  constant  equal  to  2/3, 
where  k  is  the  heat  conduction  coefficient.  The  temperature  exponent  s  is  given  by  s  =  l/2-|-2/(v— 1),  where 
V  is  the  power  index  of  the  inter-molecular  force  law.  For  argon  gas  at  STP,  v=  7.5  is  cited  by  Chapman 
and  Cowling  [6]  based  on  early  viscosity  data.  Recent  work  by  Lumpkin  and  Chapman  [7]  suggests  that 
v=  9  is  a  better  approximation,  which  is  confirmed  through  systematic  calculation  of  shock  wave  profiles. 
In  our  calculation,  the  local  Navier-Stokes  particle  collision  time  r  is  first  evaluated  according  to  r  =  /i/p, 
where  /i  ~  and  p  is  the  local  pressure.  Then,  the  new  value  r*  is  obtained  according  to  Eq.(2.7).  With 
the  general  r*,  the  BGK-NS  solver  is  used  for  the  shock  structure  solution  [12].  Since  the  BGK  scheme  is 
a  finite  volume  method,  even  with  intrinsic  unit  Prandtl  number  in  the  BGK  model,  the  heat  flux  across  a 
cell  interface  can  be  modified  to  simulate  a  gas  with  any  realistic  Prandtl  number  [12],  such  as  2/3  for  the 
current  argon  gas.  The  shock  structure  is  obtained  using  a  time  accurate  BGK-NS  solver  until  a  steady  state 
is  reached.  In  each  calculation  with  fixed  p  and  Pr,  the  mesh  size  is  chosen  to  make  sure  that  there  are  at 
least  30  mesh  points  in  the  shock  layer  and  the  whole  computational  domain  is  covered  by  200  grid  points. 


Fig.  2.1.  Comparison  of  the  theoretical  shock  thickness  \i/Ls  vs.  Mach  number  Ms  with  the  experimental  data  [2].  The 
solid  lines  are  the  results  from  the  BGK-NS  solver  [12]  and  the  new  BGK-Xu  model.  The  simulations  are  done  for  both  v=  9.0 
and  7.5  cases. 

Studies  of  the  shock  structure  are  generally  validated  by  comparing  the  reciprocal  density  thickness  with 
experimental  measurements.  The  thickness  is  defined  as 

Ls  =  (P2  -  pi)/{dp/dx)  max  • 
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Fig.  2.2.  The  density  Pn  =  (p  —  Pi)l(p2  —  pi)  distribution  vs.  xj\i  inside  the  shoek  layer  at  Ms  =  9.  Dash-dot  line, 
BGK-NS  solution  with  v=  7.5;  solid  line,  BGK-Xu  solution  with  v=  7.5;  eireles,  experimental  data  for  argon  gas  [2]. 

The  above  shock  thickness  is  normalized  by  the  upstream  mean  free  path, 

A  = 

^  Si/w  pii/2RTi 

Figure  2.1  displays  the  results,  where  “BGK-NS”  refers  to  the  solution  of  the  BGK  Navier-Stokes  solver 
with  the  original  particle  collision  time  r  =  /i/p  [12],  and  “BGK-Xu”  refers  to  the  results  from  the  same 
BGK  Navier-Stokes  solver  but  with  the  implementation  of  the  new  value  r*.  Both  v=  9  and  v=  7.5  cases  are 
tested.  All  symbols  in  Figure  2.1  are  the  experimental  data  presented  in  [2],  which  are  extensively  used  by 
many  others  to  validate  their  models  [1,  11].  The  solution  from  the  current  new  model  (BGK-Xu)  matches 
perfectly  with  the  experimental  data.  Figure  2.2  presents  the  density  distribution  p„  =  (p  —  pi)/(p2  —  pi) 
vs.  x/Xi,  where  v=  7.5  is  used  in  both  BGK-NS  and  BGK-Xu  solutions.  The  circles  in  Figure  2.2  are  the 
experimental  data  from  [2].  From  these  figures,  we  can  observe  that  the  general  particle  collision  time  used 
significantly  improves  the  results.  In  the  continuum  flow  regime,  where  the  Mach  number  of  the  shock  wave 
goes  to  1.0,  the  BGK-NS  and  BGK-Xu  solutions  converge.  In  other  words,  r*  approaches  to  r  automatically 
as  Knudsen  number  decreases. 

3.  Conclusion.  In  this  paper,  we  have  developed  a  generalized  constitutive  relation,  where  the  viscosity 
coefficient  depends  not  only  on  the  macroscopic  variables,  but  also  on  their  gradients.  Even  with  the  closure 
of  the  Chapman-Enskog  expansion  on  the  Navier-Stokes  order,  the  results  from  this  new  model  agrees  well 
with  the  experimental  data  in  the  study  of  argon  shock  structure.  The  generalization  of  the  collision  time 
from  r  to  r* , 

T 

=  l  +  T{D^glDgy 

is  important  to  capture  the  rarefied  gas  effect  in  the  continuum  transition  regime.  In  the  continuum  regime, 
such  as  the  Mach  number  approaching  to  1.0  in  the  shock  case,  the  contribution  from  (D^g/Dg)  disappears 
automatically.  This  can  be  understood  physically  that  in  the  near  equilibrium  flow  the  higher  order  con¬ 
tribution  (Burnett  D‘^g)  has  much  less  effect  than  the  lower  order  term  (Navier-Stokes  Dg).  The  further 
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application  of  this  new  BGK-Xu  model  in  the  continuum  transition  regime,  such  as  Couette  and  Poiseuille 
flows,  will  be  presented  in  subsequent  papers. 
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